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Abstract: Black holes are found to exist in gravitational theories with the presence of 
quadratic curvature terms and behave differently from the Schwarzschild solution. We 
present an exhaustive analysis for determining the quasinormal modes of a test scalar 
field propagating in a new class of black hole backgrounds in the case of pure Einstein- 
Weyl gravity. Our result shows that the field decay of quasinormal modes in such a non- 
Schwarzschild black hole behaves similarly to the Schwarzschild one, but the decay slope 
becomes much smoother due to the appearance of the Weyl tensor square in the background 
theory. We also analyze the frequencies of the quasinormal modes in order to characterize 
the properties of new back holes, and thus, if these modes can be the source of gravitational 
waves, the underlying theories may be testable in future gravitational wave experiments. 
We briefly comment on the issue of quantum (in)stability in this theory at linear order. 
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1 Introduction 

Since Einstein’s proposal in 1915, General Relativity (GR) has been established as the 
standard theory of gravitation for one hundred years. As a pillar of modern science, the 
predictions of GR have been conhrmed in all observations to date. However, one of the most 
challenging task that theoretical physicists are facing today is how GR can be reconciled 
with the laws of quantum physics to produce a consistent ultraviolet (UV) complete theory 
of quantum gravity. 

To address this issue as well as to be of phenomenological interest, extensions of Einstein 
gravity with the presence of higher order derivative terms arise in fundamental theories, 
such as string theory, loop quantum gravity, asymptotically safe gravity and others. In 
particular, it was found in [1] that adding quadratic curvature terms could improve the 
renormalizability of the underlying gravitational theory although this theory would suffer 
from an instability of ghost degrees of freedom [2]. As shown in [3], however, if the path 
integral quantization is evaluated in Euclidean space and then Wick rotated to Lorentzian 
space, this path integral can yield a theory of quantum gravity without a negative norm 
state. This approach has been applied into inflationary cosmology and provided an inter¬ 
esting interpretation for primordial perturbations [4]. 


- 1 - 




Black hole physics is believed to provide an important window to the quantum nature 
of gravity. Recently, there has been increasing interests in studying black hole solutions 
in gravity theories by taking into account higher order curvature terms, such as from the 
perspective of string theory [5-8], in the Gauss-Bonnet extended gravity [9-13], the Einstein- 
Weyl gravity [14-17], the f{R) gravity [18-21], and other cases of general quadratic gravity 
[22-26], as well as the analyses of gravitational energy of quadratic gravity [27-30]. More¬ 
over, it is important to examine the stability issue of a quadratic gravity theory by analyzing 
linear perturbations, such as in [31]. 

In addition, black hole solutions were obtained in a gravitational theory involving 
higher order terms within the scenario of asymptotical safety [32] and its stability issue 
was addressed in [33] by analyzing the so-called quasinormal modes. This study provides 
a representative example to show how the (in)stability issue of a black hole solution could 
be investigated via the method of analyzing quasinormal modes, of which the generation 
is due to the quasinormal ringing of the background spacetime under perturbations and 
hence is associated with the characteristics of black hole. The identification of the quasi¬ 
normal modes is considered to possibly falsify various black hole solutions derived in a large 
class of gravity theories through the imminent gravitational wave surveys. Note that the 
investigation of black hole perturbations has drawn a lot of interest for decades since it is 
associated with black hole stability, gravitational wave detection as well as some fundamen¬ 
tal symmetries such as the gauge/gravity duality. Analyses of these modes in theories of 
higher derivative gravity were performed in [34-37]. We refer to Refs. [38, 39] and references 
therein for recent comprehensive reviews. 

In the present paper we aim at examining the quasinormal modes seeded by a test 
scalar field propagating in black hole solutions of a pure Einstein-Weyl gravity as derived in 
ref. [16]. In Section 2 we briefly review the background gravitational theory and the associ¬ 
ated black hole solutions. Then in Section 3 we perform a detailed analysis of quasinormal 
modes seeded by a massless test scalar field that is propagating freely within various black 
hole solutions beyond Schwarzschild and investigate their behaviors. Section 4 is devoted 
to a brief discussion of the field equations for gravitational waves, in which one obtains two 
poles on the dispersion relation and hence this implies an instability of quantum fluctua¬ 
tions. We eventually summarize our results with a discussion in Section 5. Throughout the 
whole paper we use geometrized units with G = c = 1 and the (—, -|-, -|-) convention for 

the metric. 


2 Quadratic gravity and black holes beyond Schwarzschild 

We start with a brief introduction to a general theory of quadratic gravity. Consider a 
four-dimensional Einstein gravity involving quadratic curvature terms of which the action 
is written as (following ref. [16]) 

= I d\^{R- , 
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where 


C^upa — Rpi'pa 9iJ.[p^a]u 9u[p^a]p 2 9p[p9a]i/ > (^■^) 

is defined as the Weyl tensor in four dimension with R being the Ricci scalar. In the above 
action we have introduced two model parameters a and j3, which are of area dimension: 
[L^], and, describe the deviations from the regular Einstein gravity. 


2.1 Spherically symmetric and static solutions 

In order to find a black hole solution beyond Schwarzschild, we consider a static and spher¬ 
ically symmetric ansatz as follows, 

ds^ = —N{r)dt^ -\ -— -|- ) (2-3) 

F{r) 


in which two dimensionless metric factors N and F have been introduced as functions of the 
radial coordinate. As was shown in [14] as well as argued in [16], the Ricci scalar vanishes 
for any static black hole solutions of the action (2.1). As a result, the general theory of 
quadratic gravity can reduce to the pure Einstein-Weyl gravity at the classical level by 
setting (3 = 0. 

We refer to the appendix of the present paper for the details of studying a black hole 
solution beyond Einstein analytically. In the main context, we simply summarize the steps 
of constructing such a solution as follows. Firstly, one can vary the action with respect to 
the metric factors and then derive the field equations for N and F. Secondly, in order to 
exhibit the difference between this solution and the Schwarzschild one, one can parameterize 
the leading term of the F factor (denoted by Fi introduced in (A.7) in the appendix) as 
follows. 


Fi 


l±i 

rn 


(2.4) 


where rjj represents for the position of the black hole horizon and 6 is the amount of devi¬ 
ation from the Schwarzschild solution since in GR we have <5 = 0. Such a parametrization 
can provide a boundary condition for numerically solving the vacuum structure of the un¬ 
derlying gravity theory, which is the last step to implement in the whole construction. In 
the following subsection we numerically repeat the result of ref. [16] to demonstrate an 
existence of a black hole beyond Schwarzschild. 


2.2 Numerical estimates 

Note that, as has been observed in [16], for each given a, the viable value of rn is bounded. 
Specifically, if rn is too small there is no opportunity to form a non-Schwarzschild black 
hole; however, if is too large the black hole mass would become negative and hence 
leads to a quantum instability at high energy scales. For example, by setting a = 0.5, 
one numerically derives a bound: 0.876 < ru < 1.143. For any selected value of ru in 
the above bounded interval, there exists only one value of b that allows for a healthy non- 
Schwarzschild black hole. This phenomenon is related to the fine-tuning of r_|_ = 0 in the 
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Figure 1. The numerical solution of a non-Schwarzschild black hole. We choose a = 0.5, rn = 1.0, 
and finely tune <5 « 0.3633. The red, solid curve gives the numerical realization of F{r), while the 
blue, dashed curve is the numerical result of N{r). Note that, as in the treatment of ref. [16], we 
impose the normalization factor to be 0.6 for N and unity for F at infinity, in order to avoid an 
asymptotic overlap. 
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Figure 2. The comparison between the exact solution and the approximate one in (A. 11) under 
the weak field limit. The model parameters are the same as in figure 1. The upper panel shows 
the evolutions of F exactly in numerical computation and approximately in the weak field limit, 
respectively; and the lower panel describes the dynamics of N in these two cases. 


weak field limit (see the second part of the appendix). In figure 1 we show an exact example 
of the numerical construction introduced in the present subsection. 

Afterwards, using the numerical fitting one can also relate the exact numerical solution 
to the weak field limit approximate solution in (A. 10) as developed in the appendix. This 
shows that physically the weak field linearized theory can roughly describe the Einstein- 
Weyl black hole solution in the large radius regime. A specific fitting is provided in figure 2. 
One can read from the figure that the evolution of F and N in the weak field limit is in 
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agreement with the exact numerical results at large radii but deviate badly when evolving 
to the regime near the horizon. 

3 Quasinormal modes 

In this section we study quasinormal perturbations of Einstein-Weyl black holes. Quasi¬ 
normal modes have been shown to be very useful to uncover intrinsic properties of the 
geometry (for instance, see [40]). Consider a massless test scalar field ^|J propagating in 
a black hole background governed by the Einstein-Weyl theory. Its evolution obeys the 
massless Klein-Gordon equation 


O'ljj = 0 . 


(3.1) 


Plugging in the black hole metric into the above Klein-Gordon equation, one gets. 


/2F drF FdrN 


\ r ^ 2 


+ 


2N 


+ + + (3.2) 


In order to solve the above Klein-Gordon equation analytically, it is convenient to use the 
following standard separation of variables 




(3.3) 


Im 


by making use of the spherical harmonic functions. Accordingly, the Klein-Gordon equation 
can be greatly simplified as: 

2 2 rdrF rFdrN ,,, ,r‘^drF r'^FdrN-.^ ^ /o .n 

{-i:, - 1 - Yr- - + 1 ))*' + (-F- + ^ 


N 2 2N 

for each fixed value of 1. 

By introducing the generalized tortoise coordinate 


dr* = dr/\/F{r)N{r) , 


(3.5) 


we can obtain the Schrodinger-type equation for each value of I as follows, 

-Vi{r))^l{r*)=0 , (3.6) 


where the effective potential is given by 


Vi{r) = Vi{r{r*)) 


FdrN + NdrF l{l + 1)A^ 

2r 


(3.7) 


This form can be treated systematically in the analysis of quasinormal modes, which will 
be given in the following subsections. We will apply the methods of the characteristic 
integration and the WKB analysis for quantitative estimation of numerical quasinormal 
modes, respectively, in the following up subsections. 
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Figure 3. Field decay of the mode function log|<i>(t)| for a non-Schwarzschild black hole. We 
choose a = 0.7, = 1.1, and Z = 1 in numerical computations. 

3.1 Characteristic integration 

A simple but efficient way of solving 1+1 dimensional d’Alembert equations is established 
in the pioneering work [41]. In this formalism, the standard {r*,t) coordinates are replaced 
by the light-cone variables, 



u = t — r* , V = t + r* , (3.8) 

in terms of which all wave equations can have the same form. 

Considering the fact that we do not have full analytic solutions for the time-evolving 
wave equation, one efficient approach is to discretize the mode function as 

^i{N) = ^i{W) + ^i{E) - ^i{S) - ^V{S) {^i{W) + ME)) + OM) , (3.9) 

O 

where S = {u, v), W = {u + h,v), E = {u, v + h), N = {u + h,v + h), ^ and 

h is the discrete step size. Note that, in order to solve the mode function for a hxed I, 
one needs to impose the initial condition at the null boundary u = uq and v = vq. As will 
be conhrmed by the following numerical simulations, however, the characteristics of the 
associated held decay are basically insensitive to the initial conditions imposed. 

From now on we would like to drop the subscript I from the mode function but specify 
its value in detailed calculations. Through the difference equation mentioned above, one 
can get a series of time domain data <h(to)) + h), + 2h), etc, for ah possible hxed 

r*, and the quasinormal vibrations can be read-off through the transformation from the 
time domain to the frequency domain. 

The numerical computation of the evolution for a mode function |4>| is shown in hgure 3 
on a logarithmic scale. In this plot one can see a representative held decay evolution in the 
time domain. We choose a = 0.7, which ensures that the conformal term contribution is 
smaller than the minimal coupling in the action, and rn = 1.1, where exists a numerical 
non-Schwarzschild solution far from the negative mass region. One can see that the scalar 
held evolution is hrstly dominated by the quasinormal vibration, and then decays with a 
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Figure 4. Field decay of log|<i>(f)| for different background parameters. In the left panel, we 
choose a = 0.5, = 1.0, with different multipole numbers 1 = 0, 1, 2, respectively. In the right 

panel, we fix ^ = 2 and choose different values for a, which are: a = 0.3 {r^ = 0.75), a = 0.5 
(rn = 1-0), a = 0.7 {rn = 1.1). Some power law tails exist in the larger time region, but they are 
not shown due to the figure size. 


power-law tail. This is a standard scenario in the time domain prohle of black holes in 
analogue with the case of the Schwarzschild solution (e.g. see the review [39]). 

One can also investigate the held decays with different background parameters as shown 
in hgure 4. The left panel of hgure 4 shows the held decays by varying the value of the 
multipole number from 0 to 2. In this panel, one can read that for larger values of the 
multipole number I, the longer the power-law tail and the larger the slope can be in a given 
time region. This behavior is also similar to the Schwarzschild case. However, from the 
right panel of the hgure, one can explicitly read that the slope of the held decay strongly 
depends on the newly introduced model parameter a. In the limit a —>• 0, it is expected 
that one can recover the held decay of a Schwarzschild black hole. Moreover, a larger value 
of a leads to a smoother slope for the held decay for a hxed value of I as shown in the lower 
panel of hgure 4. 

We recall that there exists a bound for the model parameter a that allows for a non- 
Schwarzschild black hole solution as has been addressed in section 2.2. In the theory of 
Einstein-Weyl gravity, this bound is stronger for a smaller value of a. Moreover, this bound 
interval would move along the negative direction of the real axis with a decreasing, so we 
cannot adjust rjy as a hxed parameter. However, we hnd that the inhuence of changing ru 
(in the regular bound interval) on the properties of the held decay is not so dramatic as the 
variation of a. Also, with a smaller value of a the held decay evolution has a larger slope 
with a longer tail, which is similar to the tendency when decreasing 1. The reason could be 
understood as follows. The value of a estimates the deviation from the standard Einstein 
gravity and thus from a standard Schwarzschild black hole, which is the highest symmetric 
solution (satisfying F = N and h = 0). On the other hand, the multipole number I is 
also associated with the symmetry of a dynamical system, which is similar to the standard 
Hydrogen atom problem in quantum mechanics textbooks. Thus, the effects on the held 
decay by decreasing a and by decreasing I are similar from this perspective. 
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3.2 WKB analysis 


In order to understand the quasinormal perturbations (semi-)analytically, it is useful to 
perform the WKB analysis of the mode function. Imposing s = io; on the Schrodinger-like 
equation (3.6), one gets 

^ - (s^ + V{r*m = 0 . 

This is the Laplace transform of the original time-depended wave equation for T; 
the boundary condition 

lim = 1 , (3-11) 

r*^±oo 

which is very common in usual cases, one can get a discrete set of possible values for s(a;). 

The WKB semi-analytic approach [42, 43] is a very successful and efficient method to 
calculate quasinormal frequencies. In the present study we expand the computation up to 
the third order, and correspondingly, the square of those frequencies can be written as 

ujI = (Ko + P)- + Q) , (3.12) 

where 


and 


with ^ = n-|-1/2. In addition, the superscript (i) denotes the Tth order differentiation with 
respect to r* of the potential V{r{r*)). The subscript 0 means that the potential and its 
derivatives are calculated at the point Tq, where V{r{r*)) is an extremum. The solution of 
Lj can be determined when 

n = 0,1, 2 ... if Re(wn) > 0 , 

n =-1,-2,-3... ifRe(a;n)<0. (3.15) 

Thus we can use the formula (3.12)-(3.14) to find the quasinormal modes of the non- 
Schwarzschild black holes. 

The accuracy of the WKB method is sensitive to specific black hole solutions [44]. 
However, as has been observed in [36, 45], for low overtones (/ > n) the accuracy becomes 
better. In order to compare explicitly the method of the characteristic integration and the 
WKB analysis, we perform detailed analyses of two methods by varying the value of a from 


.(4) 


(3)' 


R = ^-'^(^ + C^)-— [ ^0 1 (7-h60^^ 

8y(2)^4+^^ 288 U(2) I + 


(3.13) 


5 Vr 


(3)4 


0 


6912 w 


(2)4 


(77 -F 188^2) - 


384 


t/(3)2t r(4) 

° ° (51 -F 100^2) -F 


1 Kn 


(4)2 


Kn 


(2)3 


2304 y(2)2 


(67-h68C^) 


+ 


1 Ko K 


(3)ta(5) 




288 y(2)2 


(19 -F 28f) - ^ 


( 6 ) 


288 w(2) 
•^0 


(5 + 4^2 


(3.14) 


(3.10) 

Under 
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a 

ru 

1 

WKB:Re(a;) 

WKB:Im(a;) 

Int.:Re(a;) 

Int.:Im(cj) 

0.1 

0.40 

0 

0.96413 

-0.98076 

1.322 

-0.784 

0.1 

0.40 

1 

1.57122 

-0.71117 

1.329 

-0.785 

0.1 

0.40 

2 

2.50942 

-0.57118 

1.901 

-0.532 

0.2 

0.60 

0 

1.25506 

-1.36612 

0.943 

-1.003 

0.2 

0.60 

1 

1.14856 

-0.58045 

1.072 

-0.567 

0.2 

0.60 

2 

1.81116 

-0.41710 

1.505 

-0.372 

0.3 

0.75 

0 

0.70972 

-0.76229 

0.848 

-0.911 

0.3 

0.75 

1 

0.94765 

-0.49453 

0.941 

-0.521 

0.3 

0.75 

2 

1.49154 

-0.37587 

1.336 

-0.343 

0.3 

0.80 

0 

0.83671 

-0.91731 

0.865 

-0.938 

0.3 

0.80 

1 

0.97668 

-0.55954 

0.950 

-0.543 

0.3 

0.80 

2 

1.51885 

-0.42566 

1.309 

-0.406 

0.4 

0.85 

0 

0.77774 

-0.84128 

0.640 

-0.543 

0.4 

0.85 

1 

0.77737 

-0.30450 

0.870 

-0.343 

0.4 

0.85 

2 

1.27841 

-0.26186 

1.209 

-0.199 

0.5 

1.00 

0 

0.54346 

-0.58875 

0.628 

-0.568 

0.5 

1.00 

1 

0.70183 

-0.29062 

0.810 

-0.369 

0.5 

1.00 

2 

1.16591 

-0.31743 

1.123 

-0.227 

0.5 

1.10 

0 

0.70962 

-0.78766 

0.705 

-0.835 

0.5 

1.10 

1 

0.81474 

-0.55477 

0.782 

-0.559 

0.5 

1.10 

2 

1.19632 

-0.37069 

1.112 

-0.399 

0.6 

1.05 

0 

0.50857 

-0.54633 

0.562 

-0.564 

0.6 

1.05 

1 

0.65280 

-0.30540 

0.768 

-0.348 

0.6 

1.05 

2 

1.05104 

-0.25422 

1.097 

-0.221 

0.7 

1.10 

0 

0.37600 

-0.39482 

0.504 

-0.495 

0.7 

1.10 

1 

0.59613 

-0.26813 

0.739 

-0.295 

0.7 

1.10 

2 

0.96098 

-0.22188 

0.710 

-0.185 

0.7 

1.20 

0 

0.47102 

-0.51351 

0.558 

-0.501 

0.7 

1.20 

1 

0.63437 

-0.36070 

0.732 

-0.333 

0.7 

1.20 

2 

0.97969 

-0.21523 

0.985 

-0.237 

0.8 

1.30 

0 

0.41946 

-0.45870 

0.539 

-0.462 

0.8 

1.30 

1 

0.54557 

-0.19984 

0.703 

-0.310 

0.8 

1.30 

2 

0.92753 

-0.25469 

0.971 

-0.187 

1.0 

1.50 

0 

0.29276 

-0.32336 

0.508 

-0.501 

1.0 

1.50 

1 

0.54461 

-0.32918 

0.653 

-0.344 

1.0 

1.50 

2 

0.83938 

-0.24689 

0.915 

-0.225 


Table 1. Values for the quasinormal frequencies for the mode function propagating in the non- 
Schwarzschild geometry based on the WKB method and the algorithm of characteristic integration, 
respectively. The value of model parameter a varies from 0.1 to 1.0 smoothly. The integer value of 
the multipole number varies from ? = 0 to 2, respectively. 
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0.1 to 1.0 and the multipole number from Z = 0 to 2, respectively. Our results are presented 
in table 1. 

From this table, we find that the results of the WKB analysis and the characteristic 
integration results are fairly consistent but not with high accuracy. Regarding this issue, we 
argue that there exists a limit on the accuracy of both two methods. For the method of the 
characteristic integration, the error by solving the differential equations could lead the wave 
function to deviate from the original form of the differential equations and the corresponding 
boundary conditions could be affected as well. To overcome this numerical deviation, one 
needs to improve the method of computer algorithm, which is a very detailed technical 
issue. Moreover, the accuracy of the results obtained from the WKB approach mainly rely 
on to which order one truncates the computation. Consider that the WKB calculations at 
higher order would become extremely lengthy and do not change the results qualitatively, 
we would like to simply take the third order truncation in our detailed analysis. 

Some generic features can be concluded here. First, in the regular bounded interval 
of all quasinormal mode frequencies have a negative imaginary part, which shows that 
scalar perturbations in the non-Schwarzschild black hole backgrounds is stable in the time 
evolution^. Second, it is observed that with increasing a, I and decreasing rn, the absolute 
value of the imaginary part of the frequency decreases. This observation is consistent with 
our previous argument on the tendency of the slope because the imaginary part of the 
frequency stands for the slope of the logarithmic time domain decay. Third, the real part 
of the frequency stands for the trigonometric vibration of the scalar field. This real part 
dramatically increases along with a larger value of the multipole number 1. 

4 Field equations for gravitational waves 

In this section we briefly discuss tensor perturbations around a non-Schwarzschild black hole 
solution. Here we use Greek letters ... to denote the coordinates on the four-dimensional 
spacetime, Latin letters i,j, k... to denote the coordinates on the two-dimensional space sub¬ 
manifold S'^ (namely, coordinates), and r, t to denote the radial and time coordinate 

respectively. And also, we use a comma to denote an ordinary derivative, while a semicolon 
denotes a covariant derivative. 

Let us consider linear perturbations about the background metric. 


t gfj,i/ -|- h 


1 


(4.1) 


where the transverse and traceless tensor satisfies 

= 0 , = 0 . (4.2) 

The first order perturbation theory will provide simple results for the perturbation of 
the Riemann tensor, Ricci tensor and Einstein tensor for spherically symmetric and static 

^We only study the parameter region of physical interest, namely, we have required a < 1 to ensure 
that the higher order derivative terms are perturbative to the Einstein-Hilbert sector, as well as a positively 
definite black hole mass. 
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backgrounds [46]. We simply summarize them as follows. By writing the components of 
the Ricci and Riemann tensors as 


Ri = g6i, R’lj=M{d’^S\-S^5i) , 

I^ = T5l, lQ=V5i, (4.3) 


where 

N{-2 + 2F+ rFr)+ rFN^r 
y - 2r‘^N ’ 

1 ^ ^ Fr 

r2 ’ ' 2riV ’ 2r ’ 

we can derive the perturbations of these geometric tensors 

6iv;, = l{-ahi + 2Mhi), 

- ViV^h{ + ViV^h\) , 

jr ;;) = ^{ViV^hj + ViV^hi - ViV^hi - ViV^h\ - vhi). 


(4.4) 


(4.5) 


Here the LHS has indexes (^, zz), while the RHS may only have {i,j). This convention 
means that the related tensors are nontrivial only when (/U, zz) are on the submanifold 8“^. 
As a result, one gets 6R = 0. 

After a very lengthy computation, one can derive the field equation for tensor fluctua¬ 
tions at leading order, which is given by 


“ + g + Rr — ^0 + + 4A4 H — )F\hi 


4 


a 


M + 


2nM^h{ 


,1 


rl 


F 


+ + + --{g - Rl) drh{ = 0 , 


(4.6) 


where we have introduced the Schouten tensor 

Pfiu = -{Rfiu — -Rdfiv) ■ (4-7) 

It is interesting to notice that, there exist the operator which appears in the first 
term of the above field equation and cannot be canceled by any constraint equations. This 
implies that at high energy the tensor fluctuations would obtain an extra degree of freedom. 
Such a new degree of freedom often corresponds to a ghost mode that would spoil the 
stability of the vacuum state quantum mechanically, such as in the Lee-Wick theory of 
particle physics. 

In this regard, the Einstein-Weyl theory of gravitation may still suffer from the quan¬ 
tum instability issue even though this theory is classically stable since the classical scalar 
perturbations can well behave as analyzed in the previous section. This instability was 
recently also addressed in ref. [47], where the authors applied the Stiickelberg approach to 
show that the interplay between the ghost graviton and the healthy graviton allows the 
theory to evade the usual strong coupling issue widely existing in massive gravity theories 
and become renormalizable, at the expense of stability. 
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5 Conclusion 


Recently, the theory of quadratic gravity has drawn the interest of theoretical physicists 
in the literature from various aspects [48-51]. In particular, it was found in [16] that 
the Einstein-Weyl gravity can allow for a black hole solution beyond Schwarzschild. In the 
present work we have revisited this type of new black hole solutions at both the background 
and perturbative levels. At the background level, we analyzed the solutions of the metric 
factors in the weak held limit and compared them with the exact numerical results. At the 
perturbative level, we have studied in detail the propagation of quasinormal perturbations 
seeded by a test scalar held within such a black hole background. 

Specihcally, we have analyzed the frequencies and time domain evolutions of quasinor¬ 
mal modes seeded by this massless scalar held in the exterior of such a non-Schwarzschild 
black hole as derived in the Einstein-Weyl theory of gravity. Our results show that the 
time domain evolution of the quasinormal modes is similar to that obtained in the regular 
Schwarzschild case where the mode functions decay in a power-law form. However, due to 
the existence of the higher derivative term, the slope of the held decay is generally smoother 
than that in the Schwarzschild one. In addition, we present a brief analysis of tensor pertur¬ 
bations, which are regarded as gravitational waves. We show explicitly that the linearized 
held equation for gravitational waves involves higher derivative operators that would bring 
the theory to be unstable quantum mechanically at high energy scales. However, it would 
be interesting to study in more depth this issue under non-perturbative approaches in order 
to reveal the relation between a ghost graviton mode and quantum renormalizability of 
gravity theories. 

A Appendix 

In the hrst part of this Appendix, we provide a detailed instruction to the background 
theory of quadratic gravity. In the second part, we show how a black hole solution beyond 
Schwarzschild can be obtained in this theory. 

A.l The theory of quadratic gravity 

Varying the action (2.1) with respect to the metric yields the background held equation as 
follows. 


Rfiu 2 -|- — 0 , 


(A.l) 


where we have introduced 



(A.2) 


which is from the conformal gravity part, and also, the Bach tensor 




fipua 


(A.3) 


due to the Weyl-Eddington term. 
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In order to study the vacuum structure of a static spacetime satisfying spherical sym¬ 
metry, we assume the absence of matter fields in the above system. As was shown in [14] 
as well as argued in [16], the Ricci scalar vanishes for any static black hole solutions of 
the action (2.1). That is, i? = 0 in the above classical theory. It turns out that, at the 
classical level, one can greatly simplify the quadratic gravity action (2.1) by taking /3 = 0, 
and therefore, the generic action reduces to a theory of pure Einstein-Weyl gravity. 

As was pointed out in [16], however, the requirement of i? = 0 does not simply lead 
to the trivial solution of the Schwarzschild black hole. This is because, by setting R = 0 
and integrating the trace of the field equation (A.l) over the spatial region could yield a 
nontrivial and non-vanishing Ricci tensor of the four-dimensional spacetime, although the 
surface term remains zero. This is also the key reason that there might exist static and 
spherically symmetric black holes over and above the Schwarzschild one as shown in [16]. 


A.2 Black holes beyond Schwarzschild 

Plugging the spherically symmetric and static ansatz (2.3) into the background field equa¬ 
tion (A.l) leads to the following two second order differential equations: 


R,rr — 


3F2 


-N - FrN^r - 


(3F -k rFr) 


2N 


Ni 


{rN^r - 2N) 

rF n 2N . .,2 Frs 1 / A^ n 

+ -r(' - + T) - a (j? - " - ’'"4 


+ 


2iv2‘'-'’ ' r ^ '’^r ' F 
r2{l-F) 2Fr.^, ,2 Fr.^, 1 

1 ^2p J 2F^ ’ 2h ’ 


rF 


2F^ 


(A.4) 

(A.5) 


where the subscript represents a derivative with respect to r. To fully determine the 
solutions of these two differential equations, one also needs to impose the horizon condition 
that 


F{r = th) = N{r = th) = 0 , (A.6) 

with rn being defined as the position of the black hole horizon. 

A.2.1 Near horizon limit 

Since the above two metric factors vanish at the horizon, in order to grasp the physics of the 
black hole solution near the horizon, it is convenient to make Taylor expansions as follows, 

OO OO 

F{r) = Y,Fi{r-rH)\ N{r) = Y,N,{r - mf ■ (A.7) 

i=l i=l 

Note that, among all the coefficients of the Taylor expansions, there exists at least one 
parameter that is a normalization factor and accordingly can be absorbed by a time re¬ 
scaling. In the following we take Ni to be the normalization factor without loss of generality. 
One can plug these expansions into the field equations (A.4) and solve for T) and Ni order 
by order. In addition, in the parametrization (2.4) we have introduced 5 to characterize the 
difference between this solution and the Schwarzschild one. 
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A.2.2 The linearized treatment in weak field limit 

From the other side, it is well known that the metric factors should approach unity when 
far away from the black hole in order to be consistent with the boundary condition of 
Minkowski spacetime. Therefore, one can analyze these metric factors in the weak field 
limit at large scales. This method has been widely applied in the literature and turned out 
to be very useful in analyzing black hole systems in modified gravity theories, for instance, 
in massive gravity models [52, 53]. In this limit, we can expand the metric factors around 
a Minkowski background as 

F{r) = 1 + /(r) , N{r) = 1 + n{r) , (A.8) 


in the limit where r is much larger than the horizon. 

Keeping leading order terms in n and /, the field equations can then be greatly sim¬ 
plified, of which the forms are given by, 

“ 4a/(r) -|- r'^n^r + ‘^r^ctf,rr = 0 , 

+ 2rn^r + r'^n^rr = 0 . (A.9) 


Consequently, the metric factors in the weak field limit can be approximately solved as 


n(r) 

f{r) 


T 

-y+T_ — 


a+mr 


+ r 


bo p (1 + mr) e 
r 2r 


r 

-mr ^ ^ 


2r 


(A.IO) 


with m? = l/2a being introduced. Moreover, the coefficients Fq, F± are integral constants 
that can be determined by numerically evolving the metric factors from the horizon to large 
length scales. 

From the approximate solution in the weak field limit, one can immediately notice that 
the appearance of the e^'^/r term would severely spoil the classical stability of the black 
hole system governed by the Einstein-Weyl theory. One possible way of circumventing 
this issue is to finely tune the value of <5 introduced in the parametrization (2.4) to ensure 
F_|_ = 0. Under this condition one gets 


/ \ 0 I "P ^ 

nir) = -h r_- , 

T T 

r 2r 


(A.ll) 


This limit also shows that we cannot choose a negative a for a regular black hole solution, 
in which case m will be imaginary and cause a vibration at large radii for f{r). Since for 
a specific numerical solution as will be constructed in the next section one can numerically 
fits the values of Fq and r_ by matching the linearized solution with the exact one. 
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